Benthos Indices: Results
Sarah Gaichas
2024-05-31
Introduction
Models for macrobenthos and megabenthos were run for spring and fall using all spatial and spatio-temporal random effects, and three catchability covariates: mean predator length at the location, number of predator species at the location, and bottom temperature at the location.
Methods
See workflow document for details.
Results
Counts of stations total and with prey
Spatial distribution of stations
# from each output folder in pyindex,
outdir <- here::here("pyindex")
moddirs <- list.dirs(outdir)
moddirs <- moddirs[-1]Spatial coverage for stomachs in fall, both groups (same set of predators, same stations)
mapfall <- paste0(moddirs[1], "/Data_by_year.png")
knitr::include_graphics(mapfall) 
Spatial coverage for stomachs in spring, both groups (same set of predators, same stations)
mapspring <- paste0(moddirs[2], "/Data_by_year.png")
knitr::include_graphics(mapspring) 
Model selection
See workflow document for details.
Model results
stratlook <- data.frame(Stratum = c("Stratum_1",
"Stratum_2",
"Stratum_3",
"Stratum_4",
"Stratum_5"),
Region = c("AllEPU",
"MAB",
"GB",
"GOM",
"SS"))# function to apply extracting info
getmodindex <- function(d.name){
# read settings
modpath <- stringr::str_split(d.name, "/", simplify = TRUE)
modname <- modpath[length(modpath)]
index <- read.csv(file.path(d.name, "Index.csv"))
# return model indices as a dataframe
out <- data.frame(modname = modname,
index
)
return(out)
}
modcompareindex <- purrr::map_dfr(moddirs, getmodindex)splitoutput <- modcompareindex %>%
dplyr::mutate(Season = modname |> map(str_split, pattern = "_") |> map_chr(c(1,2))) %>%
dplyr::mutate(Data = modname |> map(str_split, pattern = "_") |> map_chr(c(1,1))) %>%
dplyr::left_join(stratlook) %>%
dplyr::filter(Region %in% c("GOM", "GB", "MAB","SS", "AllEPU"))
benthosmax <- max(splitoutput$Estimate)
benthostsmean <- splitoutput %>%
dplyr::group_by(modname, Region) %>%
dplyr::mutate(fmean = mean(Estimate)) Do we have ecodata trends in any time series? Some.
Driven by early 80s spike?
filterEPUs <- c("MAB", "GB", "GOM", "SS", "AllEPU")
currentMonth <- lubridate::month(Sys.Date())
currentYear <- lubridate::year(Sys.Date())
if (currentMonth > 4) {
endShade <- currentYear
} else {
endShade <- currentYear - 1
}
shadedRegion <- c(endShade-9,endShade)
shade.alpha <- 0.3
shade.fill <- "lightgrey"
lwd <- 1
pcex <- 2
trend.alpha <- 0.5
trend.size <- 2
hline.size <- 1
line.size <- 2
hline.alpha <- 0.35
hline.lty <- "dashed"
label.size <- 5
hjust.label <- 1.5
letter_size <- 4
errorbar.width <- 0.25
x.shade.min <- shadedRegion[1]
x.shade.max <- shadedRegion[2]
setup <- list(
shade.alpha = shade.alpha,
shade.fill =shade.fill,
lwd = lwd,
pcex = pcex,
trend.alpha = trend.alpha,
trend.size = trend.size,
line.size = line.size,
hline.size = hline.size,
hline.alpha = hline.alpha,
hline.lty = hline.lty,
errorbar.width = errorbar.width,
label.size = label.size,
hjust.label = hjust.label,
letter_size = letter_size,
x.shade.min = x.shade.min,
x.shade.max = x.shade.max
)
fix<- splitoutput |>
dplyr::filter(Data %in% c("macrobenthos"),
Region %in% filterEPUs) |>
dplyr::group_by(Region, Season) |>
dplyr::summarise(max = max(Estimate))
p <- splitoutput |>
dplyr::filter(Data %in% c("macrobenthos"),
Region %in% filterEPUs) |>
dplyr::group_by(Region, Season) |>
dplyr::left_join(fix) |>
dplyr::mutate(#Value = Value/resca,
Mean = as.numeric(Estimate),
SE = Std..Error.for.Estimate,
Mean = Mean/max,
SE = SE/max,
Upper = Mean + SE,
Lower = Mean - SE) |>
ggplot2::ggplot(ggplot2::aes(x = Time, y = Mean, group = Season))+
ggplot2::annotate("rect", fill = setup$shade.fill, alpha = setup$shade.alpha,
xmin = setup$x.shade.min , xmax = setup$x.shade.max,
ymin = -Inf, ymax = Inf) +
ggplot2::geom_ribbon(ggplot2::aes(ymin = Lower, ymax = Upper, fill = Season), alpha = 0.5)+
ggplot2::geom_point()+
ggplot2::geom_line()+
ggplot2::ggtitle("Macrobenthos")+
ggplot2::ylab(expression("Relative benthos biomass"))+
ggplot2::xlab(ggplot2::element_blank())+
ggplot2::facet_wrap(Region~Season, ncol = 2)+
ecodata::geom_gls()+
ecodata::theme_ts()+
ecodata::theme_facet()+
ecodata::theme_title()
p
fix<- splitoutput |>
dplyr::filter(Data %in% c("megabenthos"),
Region %in% filterEPUs) |>
dplyr::group_by(Region, Season) |>
dplyr::summarise(max = max(Estimate))
p <- splitoutput |>
dplyr::filter(Data %in% c("megabenthos"),
Region %in% filterEPUs) |>
dplyr::group_by(Region, Season) |>
dplyr::left_join(fix) |>
dplyr::mutate(#Value = Value/resca,
Mean = as.numeric(Estimate),
SE = Std..Error.for.Estimate,
Mean = Mean/max,
SE = SE/max,
Upper = Mean + SE,
Lower = Mean - SE) |>
ggplot2::ggplot(ggplot2::aes(x = Time, y = Mean, group = Season))+
ggplot2::annotate("rect", fill = setup$shade.fill, alpha = setup$shade.alpha,
xmin = setup$x.shade.min , xmax = setup$x.shade.max,
ymin = -Inf, ymax = Inf) +
ggplot2::geom_ribbon(ggplot2::aes(ymin = Lower, ymax = Upper, fill = Season), alpha = 0.5)+
ggplot2::geom_point()+
ggplot2::geom_line()+
ggplot2::ggtitle("Megabenthos")+
ggplot2::ylab(expression("Relative benthos biomass"))+
ggplot2::xlab(ggplot2::element_blank())+
ggplot2::facet_wrap(Region~Season, ncol = 2)+
ecodata::geom_gls()+
ecodata::theme_ts()+
ecodata::theme_facet()+
ecodata::theme_title()
p
Trend comparisons by season, any difference macro and megabenthos?
ggplot(benthostsmean |> filter(Season =="fall"),
aes(x=Time, y=((Estimate-fmean)/fmean), colour=Data)) +
geom_errorbar(aes(ymin=(Estimate+Std..Error.for.Estimate - fmean)/fmean,
ymax=(Estimate-Std..Error.for.Estimate - fmean)/fmean))+
geom_point(aes(shape=Data))+
geom_line()+
#acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
facet_wrap(~Region, scales = "free_y", ncol = 1) +
#scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
ylab("Relative benthos biomass scaled to time series mean") +
ggtitle("Fall models: trend comparison") +
theme(#legend.position = c(1, 0),
#legend.justification = c(1, 0)
legend.position = "bottom",
legend.text = element_text(size=rel(0.5)),
legend.key.size = unit(0.5, 'lines'),
legend.title = element_text(size=rel(0.5)))
ggplot(benthostsmean |> filter(Season =="spring"),
aes(x=Time, y=((Estimate-fmean)/fmean), colour=Data)) +
geom_errorbar(aes(ymin=(Estimate+Std..Error.for.Estimate - fmean)/fmean,
ymax=(Estimate-Std..Error.for.Estimate - fmean)/fmean))+
geom_point(aes(shape=Data))+
geom_line()+
#acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
facet_wrap(~Region, scales = "free_y", ncol = 1) +
#scale_y_continuous(labels=function(x)round(x/foragemax, digits = 2))+
ylab("Relative benthos biomass scaled to time series mean") +
ggtitle("Spring models: trend comparison") +
theme(#legend.position = c(1, 0),
#legend.justification = c(1, 0)
legend.position = "bottom",
legend.text = element_text(size=rel(0.5)),
legend.key.size = unit(0.5, 'lines'),
legend.title = element_text(size=rel(0.5)))
Scale comparisons by season, any difference macro and megabenthos? Yes, way more macrobenthos, as expected.
ggplot(splitoutput |> filter(Season =="fall"),
aes(x=Time, y=Estimate, colour=Data)) +
#geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
geom_ribbon(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate, fill=Data), linetype = 0, alpha = 0.15)+
geom_point(aes(shape=Data))+
geom_line()+
#acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
facet_wrap(~Region, scales = "free_y", ncol = 1) +
scale_y_continuous(labels=function(x)round(x/benthosmax, digits = 2))+
ylab("Relative benthos biomass scaled to maximum") +
ggtitle("Fall models: scale comparison")+
theme(#legend.position = c(1, 0),
#legend.justification = c(1, 0)
legend.position = "bottom",
legend.text = element_text(size=rel(0.5)),
legend.key.size = unit(0.5, 'lines'),
legend.title = element_text(size=rel(0.5)))
ggplot(splitoutput |> filter(Season =="spring"),
aes(x=Time, y=Estimate, colour=Data)) +
#geom_errorbar(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate))+
geom_ribbon(aes(ymin=Estimate+Std..Error.for.Estimate, ymax=Estimate-Std..Error.for.Estimate, fill=Data), linetype = 0, alpha = 0.15)+
geom_point(aes(shape=Data))+
geom_line()+
#acet_wrap(~fct_relevel(Type, "Ecoregion", "Bluefish"), scales = "free_y") +
facet_wrap(~Region, scales = "free_y", ncol = 1) +
scale_y_continuous(labels=function(x)round(x/benthosmax, digits = 2))+
ylab("Relative benthos biomass scaled to maximum") +
ggtitle("Spring models: scale comparison")+
theme(#legend.position = c(1, 0),
#legend.justification = c(1, 0)
legend.position = "bottom",
legend.text = element_text(size=rel(0.5)),
legend.key.size = unit(0.5, 'lines'),
legend.title = element_text(size=rel(0.5)))
Center of gravity exploration
Has benthos shifted along the coast?
library(FishStatsUtils)
getcogVAST <- function(d.name){
fit <- VAST::reload_model(readRDS(paste0(d.name,"/fit.rds")))
dir.create(paste0(d.name,"/test"))
cogout <- FishStatsUtils::plot_range_index(Sdreport = fit$parameter_estimates$SD,
Report = fit$Report,
TmbData = fit$data_list,
year_labels = as.numeric(fit$year_labels),
years_to_plot = fit$years_to_plot,
Znames = colnames(fit$data_list$Z_xm),
PlotDir = paste0(d.name,"/test")) #already have plots, will delete this directory
saveRDS(cogout, paste0(d.name,"/cogout.rds"))
unlink(paste0(d.name,"/test"), recursive=TRUE) #removes directory with unneeded plots
}
purrr::map(moddirs, getcogVAST)Center of gravity plots with ecodata trends
plotcog <- function(d.name){
cogout <- readRDS(paste0(d.name,"/cogout.rds"))
modpath <- unlist(str_split(d.name, pattern = "/"))
modname <- modpath[length(modpath)]
cogdat <- as.data.frame(cogout$COG_Table) |>
dplyr::mutate(direction = ifelse(m==1, "Eastward", "Northward"))
ggplot2::ggplot(cogdat, ggplot2::aes(x = Year, y = COG_hat)) +
ggplot2::geom_point() +
ecodata::geom_gls() +
ggplot2::labs(y = "Center of gravity, km") +
ggplot2::facet_wrap(~direction, scales = "free_y") +
ecodata::theme_facet()+
ggplot2::ggtitle(modname)
}
purrr::map(moddirs, plotcog)## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]